Efficient estimation of nearly sparse many-body quantum Hamiltonians 
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We develop an efficient and robust approach to Hamiltonian identification for multipartite quantum systems 
based on the method of compressed sensing. This work demonstrates that with only 0(s log(d)) experimental 
configurations, consisting of random local preparations and measurements, one can estimate the Hamiltonian 
of a d-dimensional system, provided that the Hamiltonian is nearly s-sparse in a known basis. We numerically 
simulate the performance of this algorithm for three- and four-body interactions in spin-coupled quantum dots 
and atoms in optical lattices. Furthermore, we apply the algorithm to characterize Hamiltonian fine structure 
and unknown system-bath interactions. 
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I. INTRODUCTION 

The dynamical behavior of multipartite quantum systems is 
governed by the interactions amongst the constituent particles. 
Although, the physical or engineering considerations may 
specify some generic properties about the nature of quantum 
dynamics, the specific form and the strength of multi-particle 
interactions are typically unknown. Additionally, quantum 
systems usually have an unspecified interaction with their 
surrounding environment. In principle, one can characterize 
quantum dynamical systems via "quantum process tomogra- 
phy" (QPT) (Trill]. However, the relationship between rele- 
vant physical properties of a system to the information gath- 
ered via QPT is typically unknown. Alternatively, knowledge 
about the nature of inter- and intra- many-body interactions 
within the system and/or its environment can be constructed 
by identifying a set of (physical or effective) Hamiltonian pa- 
rameters generating the dynamics ll9T [l8ll . Currently, a scal- 
able approach for efficient estimation of a full set of Hamilto- 
nian parameters does not exist. 

The dynamics of a quantum system can be estimated by 
observing the evolution of some suitable test states. This can 
be achieved by a complete set of experimental configurations 
consisting of appropriate input states and observables mea- 
sured at given time intervals. Knowledge about the dynam- 
ics may then be reconstructed via inversion of the laboratory 
data by fitting a set of dynamical variables to the desired accu- 
racy. Estimating Hamiltonian parameters from such a proce- 
dure faces three major problems: (1) The number of required 
physical resources grows exponentially with the degrees of 
freedom of the system (H-d]. (2) There are inevitable statis- 
tical errors associated with the inversion of experimental data 
lUH]. (3) The inversion generally involves solving a set of 
nonlinear and non-convex equations, since the propagator is 
a nonlinear function of Hamiltonian parameters [9-18]. The 
first two problems are always present with any form of quan- 
tum tomography, but the last problem is specific to the task 
of Hamiltonian identification as we wish to reconstruct the 
generators of the dynamics. Many quantum systems involve 
two-body local interactions, so the goal is often to estimate 
sparse Hamiltonians with effectively a polynomial number of 



unknown parameters. Unfortunately, quantum state and pro- 
cess tomography cannot readily exploit this potentially useful 
feature. 

The highly nonlinear feature in the required inversion of 
laboratory data was studied in Ref.[9] in which closed-loop 
learning control strategies were used for the Hamiltonian 
identification. In that approach one estimates the unknown 
Hamiltonian parameters by tailoring shaped laser pulses to 
enhance the quality of the inversion. Identification of time- 
independent (or piece-wise constant) Hamiltonians have been 
studied for single-qubit and two-qubit cases II13[ Il4tl to verify 
the performance of quantum gates. Estimation of these Hamil- 
tonians is typically achieved via monitoring the expectation 
values of some observable, e.g. concurrence, which are time 
periodic functions. Through Fourier transform of this signal 
the identification task is reduced to finding the relative loca- 
tion of the peaks and heights of the Fourier spectrum lfl3i[l4ll . 
Bayesian analysis is another method proposed for robust es- 
timation of a two-qubit Hamiltonian II 1511 . The difficulty with 
these methods is then scalability with the size of the system. 
A symmetrization method for efficient estimation of the mag- 
nitude of effective two-body error generators in a quantum 
computer was studied in IU6I1 by monitoring quantum gate av- 
erage fidelity decay. Recently, it was demonstrated that direct 
or selective QPT schemes could be used for efficient identifi- 
cation of short- time behavior of sparse Hamiltonians 111711 as- 
suming controllable two-body quantum correlations with aux- 
iliary systems and the exact knowledge of the sparsity pattern. 
Another scheme for the determination of the coupling param- 
eters in a chain of interacting spins with restricted controlla- 
bility was introduced in Ref. IU8I1 . 

In this work, inspired by recent advances in classical signal 
processing known as compressed sensing [ 19], we use random 
local input states and measurement observables for efficient 
Hamiltonian identification. We show how the difficulties with 
the nonlinearity of the equations can be avoided by either a 
short time or a perturbative treatment of the dynamics. We 
demonstrate that randomization of the measurement observ- 
ables enables compressing the extracted Hamiltonian infor- 
mation into a exponentially smaller set of outcomes. This is 
accomplished by a generalization of compressed sensing to 
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utilize random matrices with correlated elements. This ap- 
proach is applicable for Hamiltonians that are nearly sparse 
in a known basis with an arbitrary unknown sparsity pattern 
of parameters. The laboratory data can then be inverted by 
solving a convex optimization problem. This algorithm is 
highly tolerant to noise and experimental imperfections. The 
power of this procedure is illustrated by simulating three- and 
four-body Hamiltonians for neutral atoms in an optical lattice 
and spin-coupled quantum dot systems, respectively. Further- 
more, we directly apply the algorithm to estimate Hamiltonian 
fine structure and characterize unknown system-bath interac- 
tions for open quantum systems. 



II. QUANTUM DYNAMICAL EQUATIONS 

The time evolution of a quantum system in a pure state 
is governed by the Shrodinger equation, d\4>(t)) jdt = 
— iH \ip(t)). The solution of this equation for a time- 
independent Hamiltonian can be simply expressed as = 
exp(-itH) |?^(0)). In principle, the Hamiltonian of the sys- 
tem H can be estimated by preparing an appropriate set of test 
states and measuring the expectation value of a set of 

observables {Mj} after the system has evolved for a certain 
period of time. The expectation value of these observables can 
be expressed as 

Pjk = (M 3 ) 4 , k = (i, k \ e ltH M je ~ ltH |tf fc ) (1) 

Equation (fl]i implies that the experimental outcomes {pjk} 
are nonlinear functions of the Hamiltonian parameters. To 
avoid the difficulties of solving a set of coupled nonlinear 
equations we consider the short time behavior of the system. 
Monitoring the short time dynamics of the system is valid 
when the relevant time scales of the system evolution sat- 
isfy t <C K^ 1 where, for positive operator-valued measure 
(POVM) operators {Mj}, the constant K equals 2||i7|| spec . 
The general expression of K is given in appendix B, also see 
appendix A for definition of the norms. This yields the lin- 
earized form of the Eq. (fTJ 

p jk = (V> fc | Mj |V> fe ) + it (Vfcl [H, Mj] |Vfc) + 0(K 2 t 2 ) (2) 

The linear approximation contains enough information to 
fully identify the Hamiltonian and the higher order terms do 
not provide additional information. The short-time approxi- 
mation implies prior knowledge about the system dynamical 
time-scale or the order of magnitude of | \H | | spe c- This prior 
knowledge can be available from generic physical and engi- 
neering considerations. For example, in solid-state quantum 
devices the time-scale of single qubit rotations is typically on 
the order of 1-10 ns. The switching time for exchange inter- 
actions varies among different solid-state systems from lps to 
lOOps, (for more details see appendix B.) 

We expand the Hamiltonian in an orthonormal basis {r Q }, 
where Tr(rt Tp) = d5 a ^: H = ^ Q h a T a ■ Here d is the di- 
mension of the Hilbert space. In this representation the Hamil- 
tonian parameters are the coefficients h a . The expanded form 



of the above affine equation (O is 

p jk = it^2 (ip k \ [T a ,Mj] \ip k ) h a (3) 

a 

Here we introduce the experimental outcomes as pjk — Pjk — 
(ifj k \ Mj \ip k ), since (ip k \ Mj \ifj k ) is a priori known. The re- 
lation ([5]) corresponds to a single experimental configuration 
(Mj,\ip k )). For a d-dimensional system, the total number of 
Hamiltonian parameters h a is d 2 . Thus, one requires the same 
number of experimental outcomes, pjk that leads to d 2 lin- 
early independent equations. For a system of n qubits, this 
number grows exponentially with n as d = 2 2n . In order 
to devise an efficient measurement strategy we will focus on 
physically motivated nearly sparse Hamiltonians. 

A Hamiltonian H is considered to be s-sparse if it only con- 
tains s non-zero parameters {h a }. More generally, a Hamilto- 
nian H is termed nearly s-sparse, for a threshold r), if at most s 
coefficients h a (H = ^ h a T a ) have magnitude greater than 
Vhmax where h max = max(/i Q ). By definition, the sparsity 
is basis dependent. However, for local interactions, the basis 
in which the Hamiltonian is sparse is typically known from 
physical or engineering considerations. 

III. COMPRESSED HAMILTONIAN ESTIMATION 

Our algorithm is based on general methods of so-called 
compressed sensing that recently have been developed in sig- 
nal processing theory lfl9ll . Compressed sensing allows for 
condensing signals and images into a significantly smaller 
amount of data, and recovery of the signal becomes possi- 
ble from far fewer measurements than required by traditional 
methods. 

Compressed sensing has two main steps: encoding and de- 
coding. The information contained in the signal is mapped 
into a set of laboratory data with an exponentially smaller rep- 
resentation. This compression can be achieved by randomiza- 
tion of data acquisition. The actual signal can be recovered 
via an efficient algorithm based on convex optimization meth- 
ods. Compressed sensing has been applied to certain quantum 
tomography tasks. Standard compressed sensing has been di- 
rectly used for efficient pseudothermal ghost imaging Il20ll2lll . 
Recently, a quadratic reduction in the total number of mea- 
surements for quantum tomography of a low rank density ma- 
trix has been demonstrated using a compressed sensing ap- 
proach M22I1 . 

Here, we first describe how the Hamiltonian information is 
compressed into the experimental data. The output of a single 
measurement is related to the unknown signal (Hamiltonian 
parameters) through the relation Q. Suppose we try m dif- 
ferent experimental configurations (i.e., m different pairs of 
(Mj, \ipk}))- This yields a set of linear equations 

p = ®t (4) 

where $ is a m x d 2 matrix with elements <&j k ,a = 
it/y/m(ip k \ [T a ,Mj] \ip k ) (A factor 1/y/m is included for 
simplifying the proofs, appendix C). In general m has to be 
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greater than or equal to d? in order to solve Eq. A Hamil- 
tonian estimation attempt with m < d 2 seems impossible as 
we face an underdetermined system of linear equations with 
an infinite number of solutions. However, any two s-sparse 
Hamiltonians hi and h 2 still can be distinguished via a prop- 
erly designed experimental setting, if the measurement matrix 
$ preserves the distance between hi and h 2 to a good approx- 
imation: 

(l-S s )\\h 2 -hi\\l < ||3(* 2 -Mllh < (l+^llto-Mb 

(5) 

for a constant S s 6 (0, 1). A smaller 5 S ensures higher distin- 
guishability of s-sparse Hamiltonians. The inequality relation 
(0 is termed a restricted isometry property (RIP) of the matrix 
$ l23ll . We now discuss how to construct a map $ satisfying 
this inequality, and how small the value of m can be made. 

The RIP © for a matrix $ can be established by employ- 
ing the measure concentration properties of random matrices. 
In each experiment the test state and the measurement ob- 
servable can be drawn randomly from a set of configurations 
{Mj, \ipk)} realizable in the laboratory. The independent se- 
lection of \ipk) an d Mj leads to a matrix $ with independent 
rows but correlated elements &jk,a in each row. Thus the stan- 
dard results from compressed sensing theory are not applica- 
ble here (appendix C). 

In contrast, here we derive a concentration inequality for a 
matrix with independent rows and correlated columns as the 
backbone for the RIP of our quantum problem in appendix C. 
Using Hoeffding's inequality, we show that for any Hamilto- 
nian h and a random matrix $ with column only correlations, 
the random variable ||$/i|| 2 is concentrated around \ \h\\ 2 with 
a high probability, i.e. V < 8 < 1 

Prob.{|||M||^ 2 - \\h\\l\ > S\\h\\l} < 2e- m <*V+^ (6) 

for some constants Co and c\. 

Using the above inequality, now we can show how an ex- 
ponential reduction in the minimum number of the required 
configurations can be achieved for Hamiltonian estimation. 
The inequality © is defined for any h while the inequal- 
ity in the definition of RIP, Eq.!©, is for any s-sparse h. 
As shown in Ref. 024I1 . there is an inherent connection 
between these two inequalities. It is proved that any ma- 
trix $ satisfying (O has RIP with probability greater than 
1 - 2 exp(-mc (5 f + cif + s[log(d 4 /s) + log(12e/<5 f )]). 
In addition, whenever m > c 2 slog(c? 4 /s), for a sufficiently 
large constant c 2 one can find a constant C3 > such that the 
likelihood of the RIP to be satisfied converges exponentially 
fast to unity as 1 — 2 exp(— c^m). 

The set of experimental configurations defined by Eq (4), 
and the concentration properties given by Eq (5) and (6) can 
be understood as encoding the information of a sparse Hamil- 
tonian into a space with a lower dimension. Next we need 
to provide an efficient method for decoding in order to re- 
cover the original Hamiltonian. The decoder is simply the 
minimizer of the li norm of the signal h. Implementing this 
decoder is a special convex optimization problem, which can 
be solved via fast classical algorithms, yet not stricktly scal- 
able. Furthermore, the encoding/decoding scheme is robust to 



noisy data as | \p' — |j a < e where e is the noise threshold. 
Note that e includes the error of linearization (see Eq.©) that 
is 0(^/mKt 2 ). Denote ho as the true representation of the 
Hamiltonian. For a threshold 77, ho(s) is an approximation to 
ho obtained by selecting the s elements of ho as those that are 
larger than rjh max and setting the remaining elements to zero. 
Now we state our main result: 



IV. ALGORITHM EFFICIENCY 

,4 

If the measurement matrix <£> G - s drawn randomly 

from a probability distribution that satisfies the concentration 
inequality in ((5]) with S s < y2 — 1, then there exist constants 
C2, C3, d±, d 2 > such that the solution h* to the convex opti- 
mization problem, 

minimize WhW^ 

subject to \\p' - <&h\\t 2 < e, (7) 

satisfies, 

\\h* - h \\ u < -£=\\ho(s) - h \\ h + d 2 e (8) 

with probability > 1 — 2e _mC3 provided that, 

m>c 2 slog{d 4 /s), (9) 

where the performance of a li minimizer, Eq. ([8), and the 
necessary bound 6 S < \/2 — 1 are derived by Candes in Ref. 

As an example, for a system consisting of n interacting 
qubits, the exponential number of parameters describing the 
dynamics, 2 2 ™, can be estimated with a linearly growing num- 
ber of experiments m > C2s(81og(2)n — log(s)). The sec- 
ond term, d 2 e, indicates that the algorithmic performance is 
bounded by the experimental uncertainties. Consequently, for 
fully sparse Hamiltonians and e = the exact identification 
of an unknown Hamiltonian is achievable. The properties of 
the ensemble from which the states and measurement observ- 
ables are chosen would determine the parameter S s and con- 
sequently the performance of the algorithm. The linear inde- 
pendency of the $ matrix rows for a random set of local state 
preparations and observables can be guaranteed by a polyno- 
mial level of computational overhead before conducting the 
experiments. 

A certification for the nearly sparsity assumption can be ob- 
tained from Eqs.® and (O as follows: Suppose h* m is the al- 
gorithm's outcome for m configurations. The nearly sparsity 
assumption is certified on the fly during the experiment, if the 
estimation improvement 1 1 h* n+ i — h* m \ \ converges to zero for 
a polynomially large total number of configurations. 

V. PHYSICALLY NEARLY SPARSE HAMILTONIAN 

Although physical systems at the fundamental level involve 
local two-body interactions, many-body Hamiltonians often 
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describe quantum dynamics in a particular representation or 
in well defined approximate limits. The strength of the non- 
local fc-body terms typically is much smaller than the two- 
body terms with strength J and decreases with the number k. 
For a fixed sparsity threshold 77, k v is defined as the largest 
number k for which /c-body terms have strength larger than 
77 J. Then the number of the elements of a s-sparse approxi- 
mation of a n-body Hamiltonian grows linearly as 0(ng(k v )), 
where the g(k v ) is determined by the geometry of the system. 

A general class of many-body interactions arises when 
we change the basis for a bosonic or ferminoic system ex- 
pressed by a (typically local) second-quantized Hamiltonian 
to a Pauli basis, e.g., via a Jordan- Wigner transformation. 
For fermionic systems the interactions are imposed physically 
from Coulomb's force and Pauli exclusion principle. The 
second-quantized Hamiltonian for these systems can be gen- 
erally written as: 

H = ^b pq a+d q + ^ bpgrsCL+d+ards, (10) 
p.q p,q,r,s 

where the annihilation and creation operators (dj and a~j re- 
spectively) satisfy the fermionic anti-commutation relations: 
{di, dj"} = 5ij and {a,;, dj} = l26ll . For example, in chem- 
ical systems the coefficients h pq and h pqrs can be evaluated 
using the Hartree-Fock procedure for TV single-electron ba- 
sis functions. The Jordan-Wigner transformation can then be 
used to map the fermionic creation and annihilation operators 
into a representation in terms of Pauli matrices o x , o v ,o z . 
This allows for a convenient implementation on a quantum 
computer, as was demonstrated recently for the efficient sim- 
ulation of chemical energy of molecular systems fzfo . An im- 
portant example of a Coulomb based Hamiltonian is the spin- 
coupled interactions in quantum dots which has the following 
Pauli representation: 

H = KiX-^A ® °i ® °c • • • > (U) 

where A, B,C, ■ ■ ■ indicate the location of the quantum dots, 
, a l s are Pauli operators, and fr,-^-,*,... generally represents a 
many-body spin interacting term. In practice, these Hamil- 
tonians are highly sparse or almost sparse due to symmetry 
considerations associated with total angular momentum 1 28 ] . 
For example the Hamiltonian for the case of four quantum 
dots {A, B, C, D) takes the general form JH]: 

Bexchange 

A<i<j<D 

+ (<TA-<7c)(<7B-<7d) + ((TA-O-D^B-CTc)], (12) 

Another class of effective many-body interactions often 
emerge in a perturbative and/or short time expansion of dy- 
namics, such as effective three-body interactions between 
atoms in optical lattices 12911 that we study in this work. 

Next, we simulate the performance of our algorithm for es- 
timation of such sparse many-body Hamiltonians in optical 
lattices 12911 and quantum dots 112 811 . 



A. Three-body interactions in optical lattices 

An optical lattice is a periodic potential formed from in- 
terference of counterpropagating laser beams where neutral 
atoms are typically cooled and trapped one per site. Consider 
four sites in two adjacent building blocks of a triangular op- 
tical lattice filled by two species of atoms B29I1 . The interac- 
tion between atoms is facilitated by the tunneling rate J be- 
tween neighboring sites and collisional couplings U when two 
or more atoms occupy the same site. For each site an effective 
spin is defined by the presence of one type of atom as the up- 
state f and the presence of the other type as the down-state L 
Three-body interactions between atoms in a triangular optical 
lattice can be significant. The effective Hamiltonian for this 
system is studied in Ref. 112911 . The on-site collisional interac- 
tion U, and tunneling rates J = = 2 are taken to be the 
same in all sites, also U = = C/4,4. = 2.12[/^ = lOkBz. 
The effective Hamiltonian of the 4-spin system is 

Hopt-iatt = ^ bia"af +1 + b% af a" +1 af +2 (13) 

j,a=x,y,z 

where {fr?,^} Slle functions of {J,U} and their explicit 
forms are given in appendix D. The ratio 77 = \J/U\ quan- 
tifies the sparsity level. For a fixed value of U, a smaller J 
leads to weaker three-body interactions and therefore a higher 
level of sparsity. As expected, this enhances the algorithm 
performance. 

We assume that the system can be initialized in a ran- 
dom product state \ipk) — V'fe) ® ••• ® wnere |"0fc) 
are drawn from the Fubini-Study metric induced distribution. 
The required observables for the algorithm are uniformly se- 
lected from single qubit Pauli operators {of , of, of }. This 
choice of states and observables allows for S s m 0.37 < 
V2 — 1. Let us denote the extracted Hamiltonian and the 
true Hamiltonian by H* and Btrue, respectively. Here, the 
performance of the algorithm is defined by the relative error 
1 — 1 1 if* — B t rue\\fro/\\Btrue\\fro- The results for different 
number of configurations are depicted in Fig. (Q]), for various 
values of J. As evident in Fig.([TJ, performance accuracy of 
above 94% can be obtained with only 80 settings significantly 
smaller than approximately 6 x 10 4 configurations required in 
QPT 

The robustness of this scheme was also investigated for 
10% random error in simulated experimental data leading to 
about a 5% reduction in the overall performance. 



B. Four-body interactions in quantum dots 

Another important class of effective many-body Hamilto- 
nians can be obtained for electrons in quantum dots coupled 
through an isotropic (Heisenberg) or anisotropic exchange in- 
teraction. For example the Hamiltonian for the case of four 
quantum dots (A, B, C, D) takes the general form Eq. ( TT2l . 
The first term in the summation is a two-body Heisenberg ex- 
change interaction and the last three terms are four-body spin 
interactions. In certain regimes, the ratio | J'/ J\ can reach up 
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m=number of experimental configurations 



FIG. 1: The Hamiltonian estimation average performance is illus- 
trated for a system of four adjacent sites in an optical lattice for dif- 
ferent tunneling rates, J, and collisional coupling U = 10kHz. The 
error bars demonstrate the standard deviation of the performance due 
to the random and independent selection of m configurations (shown 
only for J — 5kHz). Performance accuracy of above 90% with 
only 60 settings is achievable for J = 1kHz, which is significantly 
smaller than about 6 x 10 4 required experimental configurations in 
QPT 



to 16%. The amplitude of rj = \ J'/J\ determines the sparsity 
level of the Hamiltonian. 

Here we use an efficient modification of signal recovery 
referred as "reweighted l\ -minimization" which is described 
in appendix E. The performance of this algorithm is demon- 
strated in Fig. (O that shows a significant reduction of the 
required number of settings in contrast to the standard QPT. 
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m=number of experimental configurations 



FIG. 2: Estimation of the exchange interaction Hamiltonian for four 
electrons in quantum dots. The average performance of the proce- 
dure is illustrated for different values of \ J' / J\ with 50 iterations of 
the Zi-reweighted minimization. The standard deviations are shown 
only for \J'/J\ = 0.1 It is demonstrated that only 60 different con- 
figurations are sufficient for estimating the unknown Hamiltonian 
with an accuracy above 95% for \ J'/J\ = 0.05, instead of about 
6 x 10 4 required settings via QPT. 



VI. V. CHARACTERIZATION OF HAMILTONIAN FINE 
STRUCTURES AND SYSTEM-BATH INTERACTIONS 

A. Hamiltonian fine estimation 

In many systems a primary model of the interactions is 
often known through physical and/or engineering consider- 
ations. Starting with such an initial model we seek to im- 
prove our knowledge about the Hamiltonian by random mea- 
surements. Let's assume the initial guess about the Hamil- 
tonian Hq is close to the true form H trU e that is ||A = 
Htme — Hq\ \ <C H-fftruell- Therefore for a perturbative treat- 
ment we demand t\ | A| | <C 1, which is a much weaker require- 
ment compared to t||i/ true || <C 1. We can approximate Eq. 
(1) in the paper to find 

Pjk « (V>fc|M9|V> fe ) 

+ i(i> k \[[ e tsH °Ae- tsHo ds, M°] |^) , (14) 
Jo 

where Mf = e itH ° MjC 11110 [31]. This equation is linear 
in A, consequently, in a similar fashion as above, the com- 
pressed sensing analysis can be applied for efficient estima- 
tion of the fine structure of Hamiltonians. 



B. Characterizing system-bath interactions 

The identification of a decoherence process is a vital task 
for quantum engineering. In contrast to the usual approach of 
describing dynamics of an open quantum system by a Kraus 
map or a reduce master equation, here we use a microscopic 
Hamiltonian picture to efficiently estimate the system-bath 
coupling terms generating the overall decoherence process. 
However since we consider a full dynamics of the system and 
bath, this method can be applied to a finite size environment 
such as a spin bath, or a surrogate Hamiltonian modeling of a 
infinite bath. In the latter case a harmonic bath of oscillators 
is approximated by a finite spin bath 13211 . 

Consider an open quantum system with a total Hamiltonian: 

H = H s ®Ib+Is®H b +Hsb (15) 

and 

H SB = Yl X P<9 S P ® B 1 (16) 

where Hs (Hb) denotes the system (bath) free Hamilto- 
nian and Hsb is the system-bath interaction with coupling 
strengths { \ p . q }, and a complete operator basis of the system 
and bath being {S p } and {B q }, respectively. 

We develop a formalism to estimate A p (? parameters in the 
weak system-bath coupling regime and with the sparsity as- 
sumption that a few number of X p , q have a significant value. 

The Liouvilian dynamical equation is 

j t PSB{t) = (£o + Y, K q C pq )[p SB (t)} (17) 

pq 
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where £ [.] = ~i[H s <S) Ib + Is <& Hb, ■] and C pq [] = 
—i [S p <E> B q , .] . In the regime of weak coupling to a finite bath, 
H-HsbII "C min{||-ffs||, ||1?b||}> the Liouvillan equation (T% 
can be solved perturbatively if time t satisfies <C 1. 

For an initial system density state p k , using the matrix identity 
given in Ref. Bill we find the measurement outcomes as 



where X max means largest eigenvalue. 



Pjk « tr(p k Mj) 



(18) 



+ Y. X pM\ dse^- s ^°jC pq e sL °[pk],Mj]) 
pq J o 

where Mj is a system only observable. This affine function 
between the outcomes Pj k and coupling parameters {X pq } 
is similar to Eq.(2) in the paper for Hamiltonian estimation. 
Consequently, the compressed sensing algorithm can be em- 
ployed for computing {A P9 }s. 



VII. OUTLOOK 

We have introduced an efficient and robust experimental 
procedure for the identification of nearly sparse Hamiltoni- 
ans using only separable (local) random state preparations and 
measurements. There are a number of future directions and 
open problems associated with this work. It is not known 
how the performance of the algorithm depends on the distribu- 
tion of the ensemble from which the states and measurement 
observables are drawn. Also, a general closed-loop learn- 
ing approach for updating the knowledge of sparsity basis of 
an arbitrary Hamiltonian is an interesting open problem that 
will be of importance for generic compressed system identi- 
fication. The presented method for Hamiltonian estimation 
is promising for drastic reduction in the number of experi- 
mental configurations. However the classical resources for 
post-processing is not scalable. A fully scalable Hamiltonian 
estimation method might be achievable via a hybrid of com- 
pressed sensing and DMRG (Density-Matrix Renormalization 
Group) methods [33]. A compressed tomography method can 
also be developed for nearly sparse quantum processes (3- 



\\A\\ fro = J traced A) (A3) 



Appendix B: Analysis of the short time approximation 

The short time monitoring of the system's dynamics re- 
quires a prior knowledge of the dynamical time scales. In 
the solid-state quantum devices, in particular in the context 
of quantum control and quantum information-processing, the 
time-scale of single qubit rotations is typically on the order 
of 1-10 ns. The switching time for exchange interactions 
varies among different solid-state systems. For superconduct- 
ing phase qubit the duration of a swap gate is about 10 ns 1351 . 
For electron-spin qubits in quantum dots and in donor atoms 
(Heisenberg models) II36I438I1 . and also for quantum dots in 
cavities (anisotropic exchange interactions) ll39ll the coupling 
time is between 10-100ps, while for exciton-coupled quan- 
tum dots (XY model) and Forster energy transfer in multichro- 
mophoric complexes the relevant time scale is in the order of 
lps. Next we rigorously derive bound on the evolution time t 
that guarantees the validity of the short time approximation. 

For an input state |V>fe)> the expectation value of an observ- 
able Mj is 



Pjk = (Mt)\Mj\Mt)) = (V'fcle 



iHL 



M je 



-iHt 



m (Bi) 



Considering the expansion of the propagator e lHt = I — 
itH- \t 2 H 2 + we find 

Pjk = (i^k\Mj\ijj k )+it{tp k \[H,Mj]\^ k ) 



m\H, [H, Mj]] \tj) k ) 



(B2) 



Therefore, for the linearization assumption, it is sufficient 
to have for the Z'th term 



/ times 
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Appendix A: vectors and operator norm 

In this paper we use the following different norms: 
For a vector x, 



For a matrix A, 



\jx max {AfA) 



(Al) 



(A2) 



t l min(ip k \[H, [H, [...,M 3 ]]] \ip k ) < 

3 

t l rmn]][H, [H, [...,Mj]]]\\ 

spec 

3 



(B3) 



A tighter bound can be found for operators {Mj} from a 
POVM as 

\\[H,[H,[...,M 3 ]]]\\s P ec<2 l \\H\\ l spec (B4) 



To derive this we use 



||[A,B]|| flpec < \\AB\\ 

spec 



\\BA\\ S 



< 2 



spec ||-^|| spec 



(B5) 

and 11^11^ = 11^11^. 

This gives a single bound sufficient for linearization: t <C 

lUrrM-l 
2 1 1 II spec 
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Appendix C: RIP from a concentration inequality 

In this work, we generalize the standard compressed sens- 
ing algorithm such that the necessity for independent random- 
ness in all elements of the measurement matrix, (j>, can be 
avoided. A common approach to establish RIP ([24]) for a 
matrix $ is by introducing randomness in the elements of this 
matrix. This approach benefits from measure concentration 
properties of random matrices. In classical signal processing 
each element &jk,a can be independently selected from a ran- 
dom distribution such as Gaussian or Bernoulli. Whereas in 
the Hamiltonian estimation formulation (Eq. (4) in the paper) 
there is no freedom for independent selection of the $ matrix 
elements. 

Here we prove the concentration inequality that we em- 
ployed for establishing the restricted isometry property. 

Though $ is a random matrix, because it is constructed 
from quantum states and observables of a finite dimensional 
system, it is bounded. Thus we are able to apply Hoeffd- 
ing's concentration inequality: If vx, ...,v m are independent 
bounded random variables such that Prob.jui G [a%, &;]} = 1, 
then for S = J2i v ii 



M can be scaled by any real number, a sufficient condition is 
g/f < (1 + S)/(l — S). For the simulations in this paper, this 
ratio becomes 2.176. 



Appendix D: 4-sites optical lattice Hamiltonian 

Let us consider four sites in two adjacent building blocks 
of a triangular optical lattice filled by two species of atoms, j" 
and J.. Atoms interact by tunneling between neighboring sites, 
and J^, and through collisional couplings in the same site, 
U. The Hamiltonian for such system can be written as J29ll : 



H, 



opt—latt ~~ / 



(0.03 



Jt2 + jl2 jT3 + J i3 



u 



-0.27- 



LP 



2.l(j^ + j^j^ y y 



u 2 



Jt3 _ Jl3. 



Prob.{5 - E(S) >t}< e - 2t2 /E 4 (bi-«<) 2 
Prob.{5-E(5) < -t] < e~ 2t2/ ^ {b '' a ^ 



-0.6 



(CI) 



(Dl) 



for any t > 0. (Here E denotes the expectation value.) Set where a j ' y ' are Pauli operators 



Vi = \4>\x\ 2 for a row 
we get Vx, 



Then with S = J2i v i = \\® x \ 



Vi = x 1 '((j)i(j>\)x € (l/m)[wi,w u ]\\x\\ 2 2 
E(5)=E||$x||? e[/ !5 ]||x||? 



(C2) 



for constants wi,w u , f, g. Note that / and g are the min and 
max singular values of E(<l>^<i>). From ( IC2I ) we find Vf + . t- > 
and Vx, 

Prob.{S - g\\x\\f 2 > t+} < Pvob.{S - E(S) > <+} 
Prob.{5-/||x||f 2 < -t_} < Prob.{S*-E(S') < 

These together with dCU and (IC2b . and the choice of t + = 

(5+ and <_ = (/- l + «)||a;||? L yields 



Prob.{||$x||f 2 - ||x||fj > S\\x\\f 2 } < 2e .> a (C3) 

with e = min{l — g,f — 1}. To ensure that t+,t- > 0, 
we need 1 — #</<<?<! + <5. Since the observable 



Appendix E: Re weighted h -minimization 

In order to simulate our alogrithm performance for estimat- 
ing the above Hamiltonian we use an iterative algorithm that 
outperforms the standard l\ norm minimization B30H . This 
procedure entails initializing a weight matrix W = l^ 2 an d 
a weight factor a > 0, and repeating the following steps until 
convergence is reached: 

1. Solve for ft,, minimize HW^H^ 
subject to - <M||/ 2 < e. 

2. Update weights 

W = diagil/iM + a), l/(\h d *\ + a)). (El) 

where h = vec(hi) is the Hamiltonian vectorized form. $ is 
the measurement matrix and p' is the experimental data with 
a noise threshold e. 
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